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With the purpose of modelling the process of mountain building, we investigate the evolution of the ridge 
produced by the convergent motion of a system consisting of two layers of liquids that differ in density and 
viscosity to simulate the crust and the upper mantle that form a lithospheric plate. We assume that the 
motion is driven by basal traction. Assuming isostasy, we derive a nonlinear differential equation for the 
evolution of the thickness of the crust. We solve this equation numerically to obtain the profile of the range. 
We find an approximate self-similar solution that describes reasonably well the process and predicts simple 
scaling laws for the height and width of the range as well as the shape of the transversal profile. We compare 
the theoretical results with the profiles of real mountain belts and find an excellent agreement. 



I. INTRODUCTION 

Mountain ranges are one of the most striking features 
of the Earth and their origin and evolution have been 
investigated for a long time. It is known that the litho- 
sphere (the outer solid layer of the Earth) is a two layer 
structure in which the crust rests on the denser upper 
mantle, being separated by the Mohorovicic discontinu- 
ity (called Moho) . The lithosphere is divided into several 
approximately rigid plates that rest on the hotter and 
more fluid asthenosphere. The relative motion of these 
plates is the cause of mountain building, because of the 
shortening and consequent thickening of the crust that 
occurs when two continental plates collide (see Fig. [T] 
for a sketch) or when an oceanic plate is subducted be- 
neath a continent. On the timescale of the orogenic pro- 
cesses the lithosphere is in local hydrostatic equilibrium 
(a condition called isostasy) that implies that the visible 
regional topography is accompanied by a corresponding 
anti-topography (called root) of the Moho. 

Clearly mountain builiding is an important problem 
that involves many disciplines and interests a broad range 
of scientists. To attempt a realistic and detailed the- 
oretical description of mountain building is an exceed- 
ingly complex task (see for example the recent review 
by Avouaoi where the field data are discussed) because 
across the lithosphere there are large variations of the 
temperature, density and rheological parameters as well 
as other properties (many of which, to compound the 
issue, are poorly known). To this should be added the 
complications due to the geometry and the time depen- 
dence of the motion of the plates. Since the pioneering 
work of England and McKenzie^ several models called 
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collectively 'thin sheet models' that treat the lithosphere 
as a thin viscous layer or layers have been developed to 
take into account in a simplified way some of the above 
mentioned features (a classification can be found in Rcfs. 
iandi). These models have been used to describe moun- 
tain building, mainly by means of extensive and detailed 
numerical simulations that deal with specific ranges. 

The basic phenomena that govern the large scale evolu- 
tion of mountain belts are the spreading flow at the depth 
of the roots together with isostasy and crustal shorten- 
ing. The profile of the ridge is determined by the dynamic 
balance between buoyancy and viscous forces. Based on 
these ideas Gratton^ used dimensional arguments to de- 
rive scaling laws for the evolution of the height and the 
width of a mountain belt and argued that the evolution 
of the profile of a range is self-similar, even if he could 
not compute the exact shape. To this purpose he esti- 
mated the viscous forces assuming that the vertical gra- 
dient of the horizontal velocity takes place near the root. 
However, we shall show later that this assumption is not 
correct, since the whole lithospheric mantle is involved in 
the flow. As a consequence the scaling laws of Ref. la can 
not describe the evolution of mountain ranges. 

More recently we investigated a related problem, 
namely the formation of a ridge by the convergent flow 
of a single liquid layer over a solid moving substrate^—, 
and found that for small time T there is a self-similar 
regime in which the height and the width of the range 
scale as T 1 / 2 regardless of the asymmetry of the flow 
and the rheology of the liquid. For large time, however, 
a different self-similar regime is achieved in which the 
height and the width follow the scaling laws obtained in 
Ref. |g. Other researchers have also investigated indepen- 
dently this problem theoretically and with a laboratory 
modeli£ as well as numericall y 11 : 12 . 

Following our previous works we here reduce the prob- 
lem to its basic essentials, taking into account the two- 
layer structure of the lithosphere but disregarding rheo- 
logical and geometrical details. For simplicity we assume 
a Newtonian rheology for the crust and for the litho- 
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FIG. 1. Geometry of the two layer model employed to describe 
the formation of a ridge. The line separating the crust and 
the upper mantle is the Moho. The dashed line represents an 
isobar. 



spheric mantle, and that the problem depends on a single 
horizontal cartesian coordinate. We also ignore erosion. 
In this way we find approximate analytic solutions, scal- 
ing laws and the asymptotic behavior of the process, thus 
achieving a deeper physical understanding of the process. 

This paper is organized as follows. In section |TT] we 
describe the assumptions and we derive the governing 
equations. In section HTTl we derive the self-similar regime 
developed in the process. In section IIVI we compare the 
self-similar theoretical profile with the topography of sev- 
eral mountain ranges. Finally in sectionlVlwe discuss our 
work, whose main conclusions are: (1) the simple two- 
layer model describes quite well the evolution of many 
mountain belts, (2) their profiles have a universal shape, 
and (3) to a good approximation the evolution is self- 
similar, with the height and width increasing as T 1 / 2 . 



II. THE TWO-LAYER MODEL 

Our aim is to describe the essentials of the mountain 
building process, using a model as simple as possible, 
in order to clarify the basic physics involved. To this 
purpose we consider a two layer liquid film as shown in 
Fig. [T]and we assume for simplicity plane symmetry. The 
upper layer (the crust) has viscosity p c , density p c and 
thickness H C (X,T). The lower one (the upper mantle) 
has viscosity /i TO , density p m and thickness H m (X,T). 
Typically for a continental plate p c m 2.7 g/cm 3 , p m 
3.2 g/cm 3 , and fi c > p m . 

Initially, both layers are uniform and H c (X,0) = C 



and H m (X,0) = M. To model the basal traction that 
is believed to drive the plate motion we assume that at 
T = the bottom of the lithosphere (Z = 0) starts mov- 
ing with a prescribed velocity Ub(X). We next assume 
isostasy, which means that for < Z < H m (see the 
dashed line in Fig. [1]) the pressure does not depend on 
X . Notice that this implies that as the thickness of the 
crust increases, part of the mass of the lithospheric man- 
tle crosses the boundary between the lithosphere and the 
asthenosphere. As a consequence the mass of the litho- 
spheric mantle is not conserved. 

To derive the governing equations we assume a slow 
viscosity-dominated flow and employ a slight generaliza- 
tion of the well-known lubrication approximation (see 
for example Refs. fl^ - fl5l ) to take into account the mo- 
tion of the bottom of the lithosphere. We neglect inertia 
and assume that the slope of the free surface is gentle, 
so that the horizontal components of the velocities of the 
fluids are much larger than the vertical ones and that 
their vertical gradients are much larger than the hori- 
zontal gradients. In this way the Stokes equation takes 
the form 



dP _ d 2 U 
for < Z < H m , and 
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for H m < Z < H m + H c . In these equations P is the 
pressure, U(X, Z, T) is the horizontal velocity and g is the 
gravity. The second equations in ([1]) and ([2]) mean that 
the pressure is hydrostatic; integrating them and using 
the isostasy condition (dP/dX = for < Z < H m ) 
we find p m dH m /dX = —p c dH c /dX. Integrating this 
equation and using the initial condition we obtain 



M 



Pn 



-(C-H c ). 



(3) 



This allows elimination of H m thus yielding an equation 
for the single dependent variable H c . 

To derive the velocity profile we assume that U(Z = 
0) = Ub, that the velocity and the shear stress are con- 
tinuous at Z = H m , and that the shear stress vanishes 
at Z = H m + H c , Then we integrate twice the first equa- 
tions in ([T]) and © with respect to Z to obtain 
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0<Z< H m 

H m < Z < Hm + H r 



(4) 



Notice that the velocity profile is linear in the lithospheric mantle and parabolic in the crust and that the average 
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shear stress in the crust is exactly half of that in the 
lithospheric mantle. This means that in most situations 
the velocity drop in the crust is a small fraction of that 
within the mantle. As we will show later these features 
of the velocity field are crucial to determine the scaling 
laws for the growth of the range. 

We define the vertically averaged velocity in the crust 

as 



V, 



1 



H m +H c 



UdZ. 



(5) 



We set Ub(X) — Uqu(X), where Uq is the maximum basal 
velocity so that u verifies |u| < 1. Next we introduce the 
following dimcnsionless quantities 



h = H c /C, v = V c /U , x = X/X , t = TU /X 
Here the horizontal scale Xq is given by 



(6) 
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Finally inserting the second of dU in ([5]) and using ([5]), 
we obtain 



v = u- 



V p m M J OX M \Zflc Pm J OX 

This equation together with the continuity equation 

d{vh) 



dh 
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govern the dimensionless thickness of the crust. The pre- 
ceding equations can be easily extended to two dimen- 
sions to deal with more general geometries. 

To describe the convergence of two plates we make the 
simplest assumption: u(x) = 1 for x < and u(x) — 
— 1 for x > 0. In this way the thickness of the crust 
starts to increase in the region of convergence. The initial 
condition is h(x, 0) = 1, and the boundary conditions are 
/i(±oo,t) = 1. At x = we impose the continuity of h 
and v. 

In general this problem must be solved numerically. In 
Fig. [5] we show some solutions. All the results shown 
in figures EH and 0] were calculated for C — 30 km, 



M = 100 km, 



2700 kg/m" 3 , p m = 3200 kg/m J and 



Mc/^m = 10- These values are representative of those 
found in the lithosphere, so that the results shown can 
be applied in general to the mountain building process. 



III. SELF-SIMILAR REGIME 

We now seek the asymptotics of the problem for small 
t. We define r = (1 — p c /Pm){h — 1), then R = Cr is the 
visible topography of the range. Since at the beginning 
of the phenomenon h — 1 <C 1, the Eqs. © and ^ can 
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FIG. 2. (Color online) Numerical solutions of ((BHSJ with 
u(x) = ±1 for x ^ and h(x,0) = 1, for t =0.25, 1.31 and 
4.00 (C = 30 km, M = 100 km, p c = 2700 kg/m 3 , p m = 3200 
kg/m 3 , pc/p m = 10). 
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FIG. 3. (Color online) Evolution of the maximum relief r m = 
r(0,t). The circles correspond to the numerical solution with 
the same parameters as in Fig. [2] the straight line is r S s(0,t). 



be linearized, and with the assumption u(x) 
x -=5 reduce to 



dr dr d 2 r 
^ = ± ^ + (1 + a) 6V' 



x^0, 



±1 for 



(10) 



where a = p m C/3p c M . With typical values for the litho- 
sphere a Ri 10 -2 , so that it can be neglected and in this 
approximation the problem depends only on the scales 
X , C and Uq. 

A solution of Eg. ([T0|) as an infinite series similar to 
that given in Ref. [Z| exists (see the Appendix „). Here we 
shall show an approximate self-similar solution r ss that 
for r<l represents the asymptotics of the full solution. 
It is given by 



1 - 



Pc 
Pit. 



where 



ip erfc(^), rp = 



2y/(l + a)t 



(11) 



(12) 



Here erfc is the complementary error function. According 
to this solution the height and the width of the ridge 
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FIG. 4. (Color online) Scaled relief of the numerical solutions 



(circles) for t — 1.4 x 10 



5.7 xW° (o), 2.3 x 1(T 3 (+) 



9.7 x 1CT 2 (□) and 6.2 x 10" 1 «>)■ The solid line is r s 



follow a simple t 1 ^ 2 scaling. We define (arbitrarily) the 
dimensionless width of the ridge as w — 2x(ip = 1) = 
4y/ (1 + a)t (so the width is the distance between the 
two points in which the height is 9% of the peak height). 
Then the height and the width W — Xqw of the ridge 
are given by 



R 




W = 




p c \ Pcyw T i/2 



(13) 



(14) 



It is interesting that W depends on C, but not on Uq. On 
the other hand R depends on Uq and is nearly indepen- 
dent on C (it depends on C only through a) . Notice also 
that the aspect ratio 9 = W/R of the ridge is constant 
and equal to 



2^p c gCM 

pmUo 



(15) 



Within the uncertainties in the parameters involved, 
these formulae give the correct order of magnitude of R 
and W for real mountain ranges. 

It is interesting to compare this approximate self- 
similar solution with the numerical solutions of the full 
nonlinear problem ([8] [9]). In Fig. [3] we show the numerical 
r(0,t) and r ss (0,t). In Fig. 2] we compare the numerical 
solutions with the solution (fTTj) . From these figures it 
can be appreciated that the self-similar solution ([TT1 [T2J) 



describes quite well the shape and the evolution of the 
ridge, even for quite large t when it might be expected to 
fail (notice that the last circle of Fig. [3] corresponds to 
h(Q,t = 4) = 3, and that h(0,t= 0.62) = 1.82 for the last 
profile in Fig. [4]) . In terms of the topography this implies 
that mountain ranges whose height does not exceed ap- 
proximately 5 km are well described by ([TTr[T2")) . We then 
conclude that the self-similar solution describes reason- 
ably well the solution of the full nonlinear problem up to 
this point. We observe that for the parameters of the nu- 
merical calculations shown in these figures h(Q,t) = 4.95 
corresponds to the root of the ridge touching the asthcno- 
sphere, after which the relief can not increase anymore. 



IV. COMPARISON WITH REAL MOUNTAIN RANGES 

It is interesting to compare the present theory with the 
real profiles of mountain ranges. However at this point 
it is convenient to point out that some mountain sys- 
tems are not linear so that they can not be described 
by the present theory. For our comparisons we have 
used the digital elevation data GTOPO30 (these data 
are available in the website of the U.S. Geological Sur- 
vey's Earth Resources Observation and Science (EROS) 
Center) to obtain locally averaged profiles of 10 approx- 
imately rectilinear segments of the Alps, Andes (2 seg- 
ments), Barisan Mountains in Sumatra, Caucasus (2 seg- 
ments), New Zealand Alps, Pyrenees and Urals (2 seg- 
ments). For each segment we have drawn 50 transversal 
profiles of 101 points each. All the 10 segments we ex- 
amined have the same "pagoda roof profile. However 4 
of them (one segment of the Andes, Caucasus and Urals, 
and the New Zealand Alps) are markedly asymmetric, 
having one side steeper than the other; in addition the 
foot of the steeper side is lower than the other. 

In Fig. [5] we show the average of the 50 profiles of a 
segment of the Pyrenees along with the best fit of these 
data to af((X — b)/c) + d, where / is given in (fl~2"j) and 
a, b, c and d are constant lengths. 

In Fig. [5] we show the theoretical profile pTJTf^ and 
the 6 more symmetric average profiles. To merge these 
profiles in a single graph we plotted — di)/a,i vs. 
(X, — bi)/ci (i = 1, . . . , 6). To obtain the constants Oj, 
bi, Ci, di we followed the same procedure as we did for 
the Pyrenees. It can be appreciated that the self-similar 
approximate solution gives an excellent fit to the actual 
shapes. 



V. DISCUSSION AND CONCLUSIONS 

As can be seen in figures O and [5] the agreement of 
the profiles of actual ranges with the self-similar shape 
is very good, even for a very ancient range as the Urals. 
However some explanations are opportune. 

The theoretical profiles are sharply peaked due to the 
discontinuity of u(x) at x = 0. It is easy to solve numer- 
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FIG. 5. (Color online) Comparison of the average topography 
(dots) of a segment of the Pyrenees (shown in the inset) with 
the theoretical profile (|11H12[) . The full line is af((X-b)/c) + 
d where a, b, c and d are constants lengths determined by 
fitting the actual topography. 



ically the problem with a continuous transition of u(x). 
We have done it assuming that u — tanh(x/u>o) where 
2wq is the width of the transition. In figure |5] we com- 
pare the numerical solution for wq = 0.4 (this value was 
chosen for better visibility) with the solution for the dis- 
continuous u case for the same time (t = 1.31). We see 
that a continuous transition leads to the same profile, ex- 
cept near the top where it is rounded. The width of this 
rounded region is always ~ 2wo, but since the width of 
the range increases as t 1 / 2 the difference between the con- 
tinuous and discontinuous cases reduces with time. We 
conclude that the self-similar solution (|ll|[T2l describes 
increasingly well the profile. 

The actual topographies shown in figures [5] and [5] arc 
the result of averaging all the transversal profiles of each 
range. All the profiles employed to prepare these figures 
have a peak, but on averaging them a rounded summit 
is obtained. Notice that the noise present in the data 
due to the local topographical accidents (that occur near 
the surface of the crust and are not a consequence of the 
average lithospheric flow we are considering) introduces 
a horizontal scale u> no i SG of a few kilometers, that sets a 
limit to the size of the features that can be compared 
with the theoretical model. Then the rounded top of 
these figures whose sizes are of the order of ui no i S c do 
not contradict the sharp theoretical profile. In addition, 
this fact suggests that the transition of the basal velocity 
occurs on a horizontal scale shorter than 2u> no isc- 

In our calculations we have assumed for simplicity a 
perfect symmetry. However it is not difficult to extend 
our model to a non-symmetric situation in which |u| as 
well as C are different in each side of the ridge. To appre- 
ciate the effects of both kinds of asymmetries we show in 
figure [7] the numerical solutions for the symmetric case 
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FIG. 6. (Color online) Comparison of the theoretical profile 
from the model with the Andes (+), Caucasus (o), Alps (*), 
Urals (x), Pyrenees (□) and Barisan Mountains (CO- 




FIG. 7. (Color online) Comparison of the solution for sym- 
metric case (full line) with those for nonsymmetric basal ve- 
locity (dashed line) and nonsymmetric thickness of the crust 
(dotted line). 



and those corresponding to a nonsymmetric basal veloc- 
ity (u — 1.9 for x < and u = —0.1 for x > 0) and to 
a nonsymmetric thickness of the crust (h(x < 0, 0) = 0.9 
and h(x > 0,0) = 1.1), for t = 1.31. The parameters 
have been chosen to ensure that in the three cases the 
added dimensionless mass is equal to 2t. We can ob- 
serve that regardless of the asymmetry the crest remains 
at x = 0. For brevity we omit more details, that will be 
published elsewhere. We believe that the non-symmetric 
segments of the Andes, Caucasus, Urals and the New 
Zealand Alps can be reproduced by adequate choices of 
the parameters. 

The present theory assumes a Newtonian rheology for 
the lithosphere, although it is believed that its behavior 
is non-Newtonian. In a recent article^ we considered the 
effect of a power-law rheology in the one layer model of 
Ref. 0- We found that in the linear regime the maximum 
height and the width of the ridge increase as i 1 / 2 regard- 
less of the rheological parameters. On the other hand 
the profile of the ridge depends on the rheology, but only 
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FIG. 8. (Color online) Comparison of a profile for a continu- 
ous basal velocity (wo = 0.4, dotted line) with the discontin- 
uous case (full line), for t = 1.31. 



weakly (see figure 4 of Ref. HI) . The two layer model used 
here can be extended to include non-Newtonian behav- 
ior but to do this exceeds the scope of the present paper. 
However based on the results of one layer model we ex- 
pect that similar results will be obtained for the two layer 
model since in the linear regime both models give analo- 
gous equations. 

We do not take into account in our model the effect 
of erosion. Several authors have considered the role of 
glacial and fluvial erosion in the orogenic process, mod- 
eling the resulting redistribution of mass at large scale as 
a diffusive process (see for example Refs. [H and [13 and 
references therein). The inclusion in our model of this 
effect would modify the coefficient of the diffusion term 
d 2 r/dx 2 in equation (10). This means that a self-similar 
solution of the same kind as (11) and (12) would result, 
but with different scales. Incidentally, this could be the 
explanation why our self-similar profile describes quite 
well all the ranges analyzed regardless of their erosion 
history. Notice also that this change should not modify 
the sharp apex of the ridge, so that a rounded summit 
will not result. We leave for future work a detailed in- 
vestigation of the effects of erosion. 

The present model can be easily generalized to in- 
clude 3-D effects replacing u by a two-dimensional vector 
u = u(x, y) and d/dx by the two-dimensional gradient 
operator V = (d/dx,d/dy). The 3-D character arises 
from the dependence of u on both cartesian coordinates. 
The resulting problem must then be solved numerically. 
The 3-D effects will be important in those parts of a 
range where the average curvature radius of the crest of 
the ridge is smaller than or of the order of its width. 
On the contrary, our results can be applied whenever the 
curvature radius is much larger than the width. 

The i 1 / 2 scaling law can be justified with a dimen- 
sional argument based on isostasy, conservation of the 
crustal mass during the shortening and the balance be- 
tween gravitational and viscous stresses, entirely analo- 
gous to that employed in Ref. @. In that paper different 
scaling laws were obtained because the viscous stress was 
incorrectly estimated, since it was not realized the key 



feature of the two layer model dynamics, namely that 
the entire lithospheric mantle is involved. 

Most of the papers about mountain building deal with 
specific ranges, chiefly the Himalaya-Tibet orogeny, that 
can not be described by the present model. It is interest- 
ing to compare the results of our two-layer model with 
those of one layer-models (see for example Medvedev 12 
and Perazzo and Gratton--), and those from the two- 
layer model of Roydei>L& The one layer models consid- 
ers a single viscous layer on a solid horizontal substrate 
with convergent motion. According to Ref. and [HI 
the height and the width of the wedge increase as t@ and 
^1-/3 respectively. In Ref. it is found that (3 decreases 
with time and that the evolution of the wedge can be 
divided into three phases. Initially, /3 = 1 so that the 
wedge grows only in height. The second phase exhibits 
an almost self-similar growth in which f3 — 1/2 so that 
the height and the width increase as t 1 / 2 . For later times 
a last phase is achieved in which f3 decreases below 0.4. 
In Ref. two self-similar regimes were found correspond- 
ing to /3 = 1/2 for short times and to (3 = 1/4 for large 
time. In our two-layer model and for realistic values of 
the thickness of the lithosphere we observe only a i 1 / 2 
self-similar regime because the root touches the astheno- 
sphere before significant departures from this regime oc- 
cur. On the other hand a /3 = 1 initial phase can be 
obtained in our two-layer model if we assume that the 
the basal velocity has a continuous transition whose hor- 
izontal extent is 2wq; this phase ends around t = — Wq 
(for brevity we omit details). Thus the one-layer and 
our two-layer models yield power-laws for the evolution 
of the height and the width which have the same expo- 
nents, notice however that the factors are quite different. 

The two-layer model of Royden™ considers only the 
crust, that is divided into an upper layer with uniform 
viscosity and a lower layer in which the viscosity de- 
creases exponentially with the depth. The basal traction 
condition is assumed to hold at the bottom of the crust. 
It is shown that two regimes can occur. In the first, the 
crustal flow is directly coupled to the underlying man- 
tle. In the second the upper crustal to mid crustal flow 
is decoupled from the underlying mantle. Which one of 
these regimes occur depends on the viscosity just above 
the Moho, which in turns depends on its depth. If no sig- 
nificant low-viscosity zone develops, crustal deformation 
is coupled to the motion of the underlying mantle, and 
a triangular mountain range develops. If a low-viscosity 
zone is initially absent but develops during crustal thick- 
ening a steep-sided flat-topped plateau ultimately forms. 
If a low viscosity zone is present in the lower crust prior 
to convergence, a wide orogen with low topographic relief 
develops. In the last two cases crustal flow is decoupled 
from the mantle except at the edges of the flat region. 
The triangular profiles that are obtained in the coupled 
mode look quite similar to those obtained here. Notice 
however that the simplicity of our model has allowed to 
obtain analytic formulae for the shape of the range and 
its scaling laws, not previously known. Furthermore, 
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according to our two-layer model the flow within the 
crust should decouple from basal traction when the root 
touches the asthenosphere, being driven only by gravity, 
possibly yielding a flat-topped profile similar to those 
discussed in Ref. 0. We have not yet investigated this 
regime. 

We conclude that the simple two layer model describes 
quite well the evolution of many mountain belts. Al- 
though the lithosphere is described by many parameters, 
to a good approximation the orogenic process involves 
only Uq, C and the combination Xq (Eq. j7|). Further- 
more as long as \i m <C fi c the viscosity of the crust is 
not relevant, since most of the vertical gradient of the 
velocity occurs in the lithospheric mantle. The evolution 
of mountain belts is to good approximation self-similar 
and in the symmetric case the profile is given by (fTTI [T2l) . 
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Appendix: Linearized series solution 

Introducing the scaled variables 

t = t/(l + a), x = x/(l + a), r = r/(l + a) 

in equation (fTU|) and following the procedure described in 
the Appendix of Ref. 0, we obtain the solution for x > 

as 



2VtH- 2 (s) 



(A.I) 

where s = (t + x)j2\ft and H q (s) denotes the Hermite 
function of order q. To obtain the solution for x < one 
must change x for —x in (| A. 1|) . 
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